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ABSTRACT 



A one dimensional radiative transfer code is developed to track the ionization and heating 
pattern around the first miniquasars and Population III stars. The code follows the evolution of 
the ionization of the species of hydrogen and helium and the intergalactic medium temperature 
profiles as a function of redshift. The radiative transfer calculations show that the ionization 
signature of the first miniquasars and stars is very similar yet the heating pattern around the 
two is very different. Furthermore, the first massive miniquasars ( ^ \Q^Mq) do produce 
large ionized bubbles around them, which can potentially be imaged directly using future radio 
telescopes. It is also shown that the ionized bubbles not only stay ionized for considerable 
time after the switching off of the source, but continue to expand for a short while due to 
secondary collisions prompted by the X-ray part of their spectra. Varying spectral shapes 
also produced sizable variations in ionized fraction and temperature profile. We also compare 
the radiative transfer results with the analytical approximation usually adopted for heating by 
miniquasars and find that, because of the inadequate treatment of the He species, the analytical 
approach leads to an underestimation of the temperature in the outer radii by a factor w 5. 
Population III stars - with masses in the range of 10 — 1000 Mq and modelled as blackbodies 
at a temperature of 50000 K - are found to be efficient in ionizing their surroundings. The 
lack of very high energy photons limits the extent of heating of these first stars and have a 
distinctly different signature from that of the miniqsos. Observational effects on the 21 cm 
brightness temperature, the thermal and kinetic Sunyaev-ZeTdovich effects, are also studied 
in the context of the upcoming radio and microwave telescopes like LOFAR and SPT. 

Key words: quasars: general - cosmology: theory - observation - diffuse radiation - radio 
lines: general. 



1 INTRODUCTION 

Comprehending the end of the so called "dark ages of the Uni- 
verse" through the process of reionization constitutes a very im- 
portant step in our understanding of the formation and build up of 
structure in the Universe. It al so provides a powerful constraint on 
various cosmological models i Tozzi et alj 20001: Bharadwai & AiH 
l2005l : IZaldarriaga. Furlanetto, & Hemguisll 120041) . Observation of 
this crucial epoch in the history of the Universe is possible through 
the redshift ed 21 cm emission from neutra l hydrogen as first rec- 
ognized by ISunvaev & Ze'ldovichI ( Il975h and further improve d 
upon and developed by many authors like 'Hogan & Rees| ( Il979l) . 
IScott & Reej ( ll99(]|) and lMadau, Meiksin . &. Rees (1997). 

One of the popul ar views on how reionization occurred 
is the following (.Furlanetto. Zaldarriaga. & Hemquist 2004) ; first 
sources of light, be it Population III (hereafter. Pop III) stars, mini- 
quasars or any other other source of radiation, created ionized bub- 
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bles. These ionized bubbles probably started to form around high 
density peaks, roughly at the same time, everywhere in the Uni- 
verse. As the sources number of increased and since the expansion 
rate of these bubbles was greater than the Universe's expansion, 
the bubbles overlapped around a certain redshift. A novel class of 
large radio telescopes telescopes like LOFAk[H MW^E 21Cm43 
and are being designed, among other things, to map the 2 1 

cm emission from neutral diffuse inter galactic medium (IGM) at 
the Epoch of Reionization (EoR) as a function of redshift, where 
they will detect the statistical signature of the high redshift Uni- 
verse as a function of z ( Zaldarriaga, Furlanetto, & Hernquist 200^; 
iFurlanetto. Zaldarriaga. & Hernguisll 120041) . The ultimate goal of 
these mega radio telescopes would lie in resolving and mapping 
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the ionized bubbles and their evolution around the first sources of 
UV radiation. 

Recent observations show that quasars with black hole 
masses as lar ge as IO' ^ A/q h ave a lready existed at redshift 
6 and higher Ipan et alj I200I l2006h suggesting the existence 
of miniquasars that harbor intermediate mass black holes with 
masses in the range of 1O^~^M0 at very high redshifts (z > 
10). A scenario in which a black hole starts accreting mass 
at very high redshift thus producing UV and X-ray radiation 
that ionizes and heats up its surrounding IGM h as been stud- 
ied by many authors feicotti & Ostriker '2004a"'^; 'Nussei"2005'; 
Furlanetto, Zaldarriaga, & Hemquist 2004; Wvithe & Loeb 2004; 
Furlanetto & Loebl2002l ; lDiikstra et al.ll200 4). Many aspects of the 



above scenario are debatable, e.g., what was the primordia l black 
hole formation mec hanism jBegelman. Volonteri. & Reed I2OO6I ; 
ISpaans&Sil3l2006l) ? and what was the evolution rate of their 
mass densities as a function of redshift which var ies dramatically 
depen d ing on the model used (c.f., Beg elman. Volonterij & Reea 
l l2006h . lRhook & Haehnel^ 12006.) and .Zaroubi et aljj2007h for dif- 
ferent scenarios)? Furthermore, the spectral indices of the power 
law spectral energy distribution of a quasar, even if marginally dif- 
ferent, can suggest very different imprints. Also the accretion rate 
and lifetimes of these primordial objects are completely unknown. 

Apart from these large potential variations in parameters 
within the scenario of miniqsos, we also have the distinctively 
different scen ario of Pop I II stars a s the primary source of 
reion i zation ( Benson et al. 20061; Whalen. Abel, & NormanI 



20041; lAbel. Wise, & Brvaj|2007l; ICiardi Ferra ra, & W hitdl2003l; 
Barkana & Loebl I2OOII ; iLoebl 1200 ^; ^Chen & Mirald a-Escud3 



20061 iBromm Kudritzki, & Loebl lzQO 1 ; Kitavama et all |2004 
Wvithe & Loebll2003 V Within the context of Pop III stars there are 



many uncertainties that need to be considered. For example, the 
unknown blackbody temperature of the star (assuming they can be 
characterized by one), its mass, life-times, clustering properties 
and redshifts at which they start to appear. 

In this paper, we try to fully explore the implications of vary- 
ing the above mentioned parameters on the final maps of the 21 cm 
signal for individual objects that are initially surrounded by neutral 
IGM. In order to achieve a full sample of this parameter space it 
becomes imperative that a fast and realistic method of modelling 
these parameters is introduced. Previous simulations of radiative 
transfer in cosmological context and studies on the impact of first 
sources on the surrounding IGM, show that the regions that are 
ionized around these first sources, ar e almost spherically symmet- 
ric, e.g.. Figure 7 of lKuhlen&lladau ( 2005) . Thus we resorted 
to developing and using a 1-D radiative transfer code to study the 
ionizing fronts (I-front), their velocity, size, shape and so on. 

Quick execution times and the simplicity in setting up the ini- 
tial conditions allow us to simulate the influence of the miniquasars 
with varying masses for different redshifts, lifetimes and spectral 
indices. Specifically, we explore the influence of individual black 
hole masses ranging from lOOM© to 10* Mq, redshifts between 
10 and 30 and lifetimes in the range of 3 to 20 million years. The 
influence of first stars on the IGM was explored with mass ranges 
of IO Mq-IO'^Mq, w hose luminosity is determined mainly by their 
mass jSchaereJl2002l) . 

The upcoming radio telescopes like LOFAR, MWA, 21CMA 
and SKA will be used to observe the redshifted 21 cm from the 
dark ages and the epoch of reionization. Hence, trying to predict 
the topology, intensity and distribution of the ionized regions as a 
function of redshift is crucial in the development of sophisticated 
data analysis techniques to retrieve the 21 cm EoR signal. Given 



the spatial resolution of these radio telescopes, it is very hard to 
differentiate b etween the ionization bubbles ar o und miniqsos and 
Pop II I stars. IC hen & Miralda-Escude (2006^; 'Kuhlen & Madau| 
(2005l); IZarou"bi & Silk (2005) have suggested ways to differentiate 
between various sources of ionization. In this paper it is shown that 
the heating caused by these two possible ionization sources is sig- 
nificantly different especially at the early stages of the EoR, leading 
to different strengths of coupling between the spin and kinetic tem- 
perature which is reflected in the measured brightness temperature. 
Therefore we argue that this feature could, in principle, be used to 
distinguish between ionization sources. 

The paper is structured as follows: In section ij2]we detail the 
rate equations used in the simulations and also describe how we 
solve for the radial dependencies of the ionized fractions of the 
various species involved and of the temperature. Accuracy and re- 
liability of the radiative transfer code developed were tested. Re- 
sults of the two tests conducted, namely, propagation of the ion- 
izing front into the mean non-expanding IGM and the Shapiro & 
Giroux case in which the IGM is expanding with the Universe, are 
presented in fj3] Sj4] outlines some of the results obtained in the ap- 
plication of this code to power-law sources. A similar exercise is 
carried out for the case of stars and results discussed in f|5] The cou- 
pling mechanisms between the spin and kinetic temperatures and its 
effect on the brightness temperature is discussed in Sj6] Other ob- 
servables like the thermal and kinetic Sunyaev-Ze'ldovich effects 
caused by these sources are considered in SjT] We present our con- 
clusions and the possibilities of expanding this work further in ijs] 
We also present a comparison of the heating profiles provided by 
the code developed here with an analytical approach developed ear- 
lier in the appendix [jA] 



2 SIMULATIONS 

The results presented in this study are based on a 1-D radiative 
transfer (RT) code with which the evolution of HI, HII, Hel, Hell, 
Helll, free electrons and t emperature is monitored in tim e. In order 
to achieve this, we follow [pukugita & Kawasaki' ('19941), wherein a 
set of coupled differential e quations linto[3landll2ll are so lved . We 
follow I Jones &Wvsel fl985 ) and Fukugita & Kawasa^ h994h in 
treating hydrogen as two-level system plus a continuum. 

For this study, in the case of miniqsos, we simply assumed that 
the density surrounding it is the mean density in an expanding Uni- 
verse. For the case of stars, we assume an isothermal like profil e of 
the for m p(r) [cmr'-^] = 3.2 X (91.5 pc/r)*^, as in lMellema et al.l 
( l2006l) . until the density reaches the mean IGM density, after which 
we set the density to that of the mean IGM. The cosmological pa- 
rameters used here are those set by the WMAP 3'^'' year results 
(iSpergel et al.2007,) . 



2.1 Rate equations 

In the RT code we developed, the radiation is allowed to propagate 
radially with an underlying density profile, which in our case is the 
constant background IGM density. The code solves the following 
rate equations and finds the number density of each species in time 
and space: 



diiH„ 
dt 



= FHinHi — ctHiincnHii 



(1) 
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where, nHj , nHu , nHcj , nnon and nnoni are the neutral & 
ionized hydrogen, neutral, single and doubly ionized helium den- 
sities respectively. Uo is the total electron density given by nc — 
nH„ + HHei + 2nHG„ • Eh J , Ehsj and Ehbu are the ionization 
energies for the corresponding species ./? are the coUisional ion- 
ization coefficients, a and ^, the recombination coefficients, a are 
the bound-free photo-ionization cross-sections, 72c ~ ctn, (T7) x 
(meA:r^/27r)^''^e~^"''^/^^ , is the photoionization coefficient due 
to background photons. The subscripts denote the species to which 

the coefficient belongs. 

The radiation flux N{E\ r; t), is the same as in lZaroubi & Silkl 

( l2005h : 



N{E-r-t)^, 



-T{E;r;t) ^9 

(r/Mpc) 



:HE) [c 
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where I{E) is the spectral energy distribution, Ag, is the normal- 
ization coefficient calculated as; 



Ag 



Etotal 



J I(E)dE' 



(8) 



where Etotai is the total energy output of the^Vit objects within the 
energy range (Erange), and the optical depth t{E; r; t) is given by; 



T{E-r-t) = ^ j a,{E)n,{r-t)dr 



(9) 



Where the sum is over all species. The effects of secondary ion- 
izations - due to the kinetic energy carried by ejected electron - 
have been folded in through the terms fn and fne- These two 
terms depend of the io nization state of the medium and given by 
jShuU & van Steenbera.1985.) ; 



fH = 0.3908(1 
and. 



.4092-. 1.7592 
•^ion } 



fHe = 0.0554(1 - X. 



.4614n1.6660 



(10) 



(11) 



Here Xion is the ionized fraction of hydrogen. For ionization, 
this effect becomes significant when dealing with very high energy 
photons ( > IKeV) and when the medium is neutral to partially 
ionized. 

Note here that the causality aspect of the code comes in 
through the ionization terms. Because, as will be discussed later, 
we are presenting a grid code in which the memory of the ioniza- 
tion history of all the cells prior to the cell under consideration is 
embedded in the optical depth, which in turn dictates the amount of 
radiation available to ionize a particular cell. 

Fits for the coUisional ionization and recombination coeffi- 
cients were obtained from the appendix of iFukugita & Kawasakil 
( Il994l) . The photoioniza tion cross-sections t hough are obtained us- 
ing the fitting formula o f lVemeretal.l ( ll996h . For the exact form of 
these fitting function we refer the reader to these two papers. 

The temperature evolution is monitored by coupling equations 
[T]to[3]with the equation below. 
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(12) 



The form of equatio n [72l is identical to that of 
iFukugita & Kawasakil ( Il994h except for the inclusion of the 
Compton heating term (ref. iMadau & Efstathi ou 1999) which 
become important as we approach the source which is placed at 
the centre, i.e., radius equals zero. In the above equation, Ei is the 
threshold energy of ionization of the i*'' species which is either 
HI, Hel or Hell, as, is the Thompson scattering cross-section 
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of an electron (cts — 6.6524 x 10~^^ cm'^). and rji are the 
coUisional-ionization cooling and recombination cooling coeffi- 
cient, respectively. LdHej, is the dielectronic recombination cooling 
coefficient due to Hen. i^i is the collisional excitation cooling 
coefficient and On the free-free cooling coefficient. Cooling due to 



Hubble expansion is accounted for by the term, ^ y — ^T^j where 
^ = // is the Hubble constant. The factor fHeat is the amount of 
heat deposited by secondary electrons and is given by; 



fHeat — 



0.9971(1 - (1 ~ xtf'^^y-'''"'') 
0.15 



if Xion 
if Xion 



> 10- 

s; 10" 



Th e above expression is an ext rapolation of fitting formula 
used in IShull & van Steenberd and the results of their cal- 

culations as plotted in their Figure 3. We here use the fact that the 
heating fraction never really goes to zero but saturates around a 
value of 0.15. The fitting formulas used in lShull & van Steenberd 
( Il985h are appropriate for relatively high energies (typically > 100 
eV). 



fHixio„,E) ^0.3908(1 
where a(xion, E) is, 



.4092a(a;io„,E)Nl.7592 



a{xion,E) = — arctan 



V0.12KeVy 



0.03 , 



(13) 



(14) 



Notice that the inclusion of an energy dependent fitting for- 
mula for the fraction of energy that goes to ionization, as in equa- 
tion [T3](Diikstra, Haiman, & Loeb 2004), does not alter the result 
significantly. This is due to the fact that the lower energy photons 
are trapped close to the vicinity of source for ionization and the re- 
maining p hotons are of relatively high en ergies for which the fitting 
formula of'Shull & van Steenberg ( 1985) holds. 

The integrals in equation [T] to [T2] are pre-calculated and tab- 
ulated as a function of UiXi, where the index 'i' refers to HI, Hel 
or HelL n; refers to the abundances and Xi the ionized fraction 
of the ith species. In other words we are tabulating the integrals 
as a function of optical depth per unit distance. Therefore in prin- 
ciple we should have a three dimensional table, as is commonly 
used jlliev et Jl2006l) . But the integration table in our implemen- 
tation is at most two dimensional. We get around using three di- 
mensional tables by utilizing the fact that the functional form for 
the photoionization cross-sections of Hi and Hen are the same ex- 
pect for a constant differen ce (see equations (B13) and (B16) of 
iFukugita & Kawasaki|[T99i) . 



2.2 The Algorithm 

The aim of the radiative transport code developed is to compute 
the fraction of the ionization for the species of hydrogen and he- 
lium and the temperature for points along a radial direction away 
from the source, at various times. This was achieved by solving the 
time dependent rate equations (eq uationsllll2t using the ODEINT 
routine of Numerical Recipes ( Press et al]| 19921) . using the implicit 
scheme of integration for stijf equations. The entire code was de- 
veloped in ANSI C. 

Figure[T]shows the basic steps involved in the solution of these 
rate equations. Firstly, the direction radial to the source is gridded 
into cells, the sizes of which are decided as described later. The 
inputs to the code at this stage are the parameter space we desire 
to probe in this study which includes the spectrum, masses of the 
sources and the redshifts at which these sources switch on. The 



^ start ^ 



Grid the radial direction from 
Rstart - Rmax with resolution Ar 



I 



/ Input 

' Type : (Star/AGN) 
Mass, Lifetime, Redsliift / 



Setup iuitiai conditions for the 
various species and temperature 
SET: time = 



Generate the loolc tables for ttie 
integrals in the rate equations 



SET Grid Index 



Solve Rate Equations for Grid 'i 
and evolve it for a period At 
THEN:i = i+I 




Figure 1. Flowchart outlining the basic modules in the algorithm imple- 
menting the radiative transfer code. 



spectrum of the sources, broadly speaking are classified into power- 
law and blackbody. The slopes of the power law and the tempera- 
ture of the blackbody can be set to any desired value. Once the 
spectrum of the source is known, a lookup table for the values of 
the integrals appearing in the rate equation is created. 

Before the start of the simulation, the initial conditions are 
setup as follows. The region around the source is set to be com- 
pletely neutral with the primordial abundance fraction of hydro- 
gen and helium. All species therefore have densities correspond- 
ing to that of the IGM at that epoch, nn{z) = nH(0)(l + z)^ , 
n-He(-z) = nHo(0)(l + 2;)^, wherenH(O) = 1.9 x 10"'^ cm"^ and 
?iHe(0) = 1.5 X 10~* cm"'' are the IGM densities at redshift zero 
of hydrogen and helium respectively or the initial density profile 
could also be set according to a profile as in equation [21] Abun- 
dances of all singly and doubly ionized species riHn , nncn ^^Hcni 
and electron density rie is set to zero. Initial kinetic temperature 
To = TcMB (2) (1 + 2)^/(1 -I- 250). Here we have assumed that the 
kinetic temperature of the gas was coupled to the CMB temperature 
(Tcmb) till a redshift of 250, after which the Tcmb continued to 
fall off as 1/(1 + z) and Te like 1/(1 + zf. 

We start our simulation at Rstart, typically 0.1 physical kpc 
from the location of the source. All hydrogen and helium are as- 
sumed to be completely ionized inside this radius Rstart - The out- 
come of the calculations for different choices of Rstart has been 
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tested and, for the choice made in this paper, found to converge to 
the required accuracy. 

Each cell is then updated for time At, which again is chosen 
based on a convergence criterion. This At is not the intrinsic time- 
step used to solve the differential equation itself because that is 
adaptive in nature and varies according to the tolerance limit set in 
the ODE solver. On the other hand, the At here decides for how 
long a particular cell should evolve before moving on to the next. 
For example, we cannot evolve the first cell for the entire lifetime 
of the source and move on to the next cell. Therefore, At is decided 
by reducing it by half each time until two consecutive runs with At 
and At/2 give the same final result for the ionized sphere within a 
tolerance limit. As eluded to before, the code is causal in the sense 
that a cell i+1 is updated after cell i. Note that the light travel time 
is not taken into consideration explicitly since the ionization front 
(I-front) is very subluminal. Hence, all cells are updated for time 
nAt at the n*'^ time-step. 

After all cells until the last cell iRmax are updated for time 
At, the resulting values are stored and then passed on as initial 
conditions for the evolution of the cell in the next At interval of 
time. The update of all Ucm is repeated n times such that nAt = 
tquasar, whcrc tquasar IS the life time of the miniquasar. The various 
quantities of interests can be stored in a file at intervals of choice. 

A radial coverage of Rmax is chosen a priori which, depend- 
ing on the problem, can be set to any value. Typically we do not 
need to go above ten comoving mega-parsecs. This radius is then 
gridded equally with a resolution of Ar, which like the time resolu- 
tion, is decreased to half its value until it meets a given convergence 
criterion, which here is that the final position of the I-front is accu- 
rate within 0.5%. For a typical run, which normally takes about 100 
time steps, this criterion gives an accuracy of about 5 x 10~^ per 
time step. 

The block computing the solution of the rate equations in Fig- 
ure [T] requires further explanation. Once the spatial resolution is 
set by following the procedure described above, we compute for 
every cell, ^{niXi)s, where j is the grid-cell under consider- 
ation, and these computed sums are used to evaluate the values of 
the integrals in equation[T]to[T2]from pre-computed tables through 
a simple polynomial interpolation. The initial conditi ons and the 
value of the integrals are passed to the solver ODEINT jPress et al] 
Il992h with the driver to solve the stiff equations. At the end of this 
run all cells are updated for a time At. 
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Figure 2. Top panel show.s the simulated (.solid) and the theoretical (dashed) 
results of the evolving I-front in an uniform background. This simulation 
include only hydrogen of constant density (i.e., no expansion). 



3.1 Test-1: I-front velocity in simple cosmology 

The first problem was to test the expanding HII bubble around a 
source that produces a fixed number of photons per unit time. Ana- 
lytically, the position ri and velocity vi of the I-front can be written 
as: 



: rs[l - exp{-t/trcc)] 



1/3 



(15) 



exp( — t/trcc) 



3t,ec [1 - exp{~t/tre.)]y3 ' 



(16) 



Where rs is the Stromgen radius jPopita & Sutherlanci2003h ; Ts 



= [CaB{T)nB] ^, is the recombina- 



3iv-, 

tion timescale; nn, the neutral hydrogen density; C is the clump- 
ing factor and afl(T) = 2.6 x 10""(r/10'*)"°-*'^ cm^s"^ is 
the hydrogen recombination coefficient at temperature T. The pa- 
rameters used for this example are: N-y = 10^* S^^, (7 = 5, 
Qs = 2.6 X 10"" cm^s-^ and nn = 1.87 x 10"" cm"^. 

The solid line in Figure|2]is plotted based on equation [T6l and 
the dashed line represents the numerical solution to the problem. 
The numerical solution obtained by the RD code is within 0.1% of 
the theoretical value. 



3 TESTING THE CODE 

In order to examine the performance of the code we carry out three 
tests. The first test is against the analytical form of the velocity 
of the I-front with constant background density and central source 
emitting a fixed number of photons, namely, the Stromgen sphere 
case. The second, is to compute the I-front velocity around a source 
again with a constant photon flux but with the background density 
evolving with redshift. This i s then compared with the analytical 
solution o f lShapiro & GirouxHl987h . In both these cases the I-front 
is defined as the position at which the ionized and neutral hydrogen 
fractions are equal. The third test was a comp arison between an 
analytical model of heating tZaroubi et al.ll2007i) and the RT code. 
The results agree well except in the discrepancy are discussed in 
the Appendix. 



3.2 Test-II: The Shapiro & Giroux test: Ionization front in 
an expanding Universe 

Shapiro & Giroux (1987) derived an analytical solution for the po- 
sition of the I-front as a function of cosmic time in an expanding 
Universe that contains only hydrogen. The solution is of the form, 

y{t) = Ae^'""[t/t.£2(At./t) - EiW], (17) 

where j/(t) = {ri{t)/rs.i)^ , ti is the age of the Universe at the 
time the source has turned on and rj and rs,i are the comov- 
ing I-front position radially from the centre of the source and the 
initial Stromgen radius {3N.y/CaB{T)n'jjy^^, respectively. Here 
uh is the comoving density of hydrogen at source turn-on time. 
A = tiCaB{T)nH, is basically the ratio between cosmic source 
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Figure 3. Top panel shows the simulated (solid) and the theoretical dashed) 
results of the evolving I-front in an expanding Universe with background 
density being that of the IGM at that epoch. 



turn-on time and the recombination timescale. £"2(2;) is an expo- 
nential integral of the form £2(2) — e'"^^ /t^dt. 

Result of the test is plotted in Figure [3] The numerical and 
analytical solutions here also agree to within a couple of tenth of a 
percent throughout the evolution of the I-front. 



4 APPLICATION TO POWER-LAW SOURCES 
(MINIQSOS) 

Sources of ionizing radiation with a power-law spectral 
energy distribution has been co nsidered by many authors 
(e.g ., Madau, Meiks in!"& ReesI 



('1994'); 



1997): Kuhlen & Mada tj feOOSh: 
Nussei< ( 120051) ; IZaroubi & Silkl 



[Fukugita & Kawasaki 
(2005); Zaroubi et al. (2007)). The radiative transport code is ap- 
plied to these power-law sources where the ionization and heating 
patterns around them are studied as a function of their spectral in- 
dex, mass and redshifts. For some of the plots in this section we 
have assumed masses of quasars in the order of IO^Mq. Although 
recent works like Volonteri, Lodato, & NataraianI j2007h does show 
that there could be relatively massive blackholes at early redshifts, 
a 10'' M0 would be unlikely. But we have included this case in our 
discussion for a completness in spanning the parametre space and 
also to clearly distinguish the various effects of these sources on 
the IGM. In the following subsections we introduce and discuss 
our main findings. 



4.1 Energy spectrum the miniquasar 

Rece nt o bservations and catalog ues published in the litera- 
ture dVanden Be rk et al. 2001; Vig nali, Brandt. & Schneideill2003l ; 
iLaor et al.H 19971 ; lElvis et al][l994l) . suggest that the energy spec- 
trum of quasars follows a power-law of the form E^" . Specifically, 
we explore two types of power law spectra: 



(18) 



The value of a is fixed to unity for most of the study, although 
in section we do briefly discuss the effects of varying this pa- 
rameter. LE and HE stands for the 'Low Energy' and the 'High 
Energy' lower limits to which the spectral energy distribution ex- 
tents. The HE case is considered in order to take into account, in 
an approximate manner, the possibility of the lower energy pho- 
tons being absorbed in the close vicinity of the source. Instead 
of single slope one c an also adopt multi-slope spectral templates 
dSazonov. Ostriker. & SunvaevI 120041) but this is not done in this 
study. 

The miniquasars are assumed to accrete at a constant fraction 
e (normally 10%) of the Eddington rate. Therefore, the miniqso 
luminosity is given by; 



L = eLodd(M) 



1.38 X 10'' 



M 



erg s 



(19) 
(20) 



1(E) oc E" 



if 10.4eV <E < IkeV (LE case); 
if 200eV <E < IkeV (HE case). 



The luminosity derived from the equation above is used to nor- 
malize the relation in [TS] according to equation [S] The normaliza- 
tion, both for the HE and LE case, is done for an energy range 
of 10.4 eV to 1 keV. Simulations were carried out for a range of 
masses between 10 and lO'' M©. Although the number of photons 
at different energies is a function of the total luminosity and spectral 
index, if we assume that all photons are at the hydrogen ionization 
threshold, then the number of ionizing photons thus obtained for the 
mass range given above is in the order of 10^" to 10^^ These are 
of the same order of magnitude of the number of io nizing photons 
being e mployed for simulations by various authors. iMellema et al.l 
( l2006l) . lKuhlen & Madaul i2005h . for example. 



4.2 Radial profile of species' fraction 

Plotted in Figures |4] and |5] are snapshots of the neutral fraction as 
a function of the radial distance away from the centre of the black- 
hole for two different cases. The first (FiglJl is the HE case, i.e., 
without low energy photons. Whereas, the second (Fig|5} is for the 
LE case , i.e., the photon energy range spans 10.4 eV to 10'* eV. 
As indicated on the Figure, each panel corresponds to a different 
mass. The lines in each panel corresponds to snapshots of the neu- 
tral fraction at 1, 3, 5 and 10 Myr. In both cases the sources are 
turned on at a redshift of 20. 

Figure[5]clearly shows that the size of the ionized regions are 
substantially larger than those in Figure|4] Even the 1000 solar mass 
miniqso produce regions that are tens of kpc in physical size. This 
is due to the increase in the number of photons in the lower energy 
(LE) range (13.6 to 200 eV), which increases the probability of a 
photon - atom interaction relative to the HE case (cr oc E^'^). Equi- 
librium is only reached when the mini-quasar life-time approaches 
the recombination time-scale. Since the problem is set in an ex- 
panding Universe we nev er actually reach an equilibrium solution 
(IShapiro & Girouxlll987h . We see that the ionized bubble gets to a 
couple of hundred kpc for the case of high energy photons and to 
more than 0.5 Mpc for ICf' Mq miniqsos with low energy photons. 
All distances plotted in the Figures are physical unless specified 
otherwise. 

The abundance of HI, HII, Hel, Hell and Helll as a function 
of radial distance are tracked in time. Figures |6] and |7] shows the 
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Figure 4. The neutral hydrogen fraction as a function of distance from the 
centre of the black hole with masses 10^'*'^'® Mq (top-left to bottom right 
panel) at a redshift 2 = 20 is shown. Four lines in each panel indicates the 
position of the I-front after 1,3,5 and lOMyr. These are miniquasars with no 
UV ionization photons. (200eV < E < lO^'eV). 



fraction of these species after 10 Myr of evolution as a function of 
radius. The miniquasar is switched on at redshift of 20. 

The interplay between the evolution of these species among 
each other and with the gas temperature provides us with rich and 
interesting structures in profiles. Miniquasars of 10 and 100 solar 
masses are not able to produce substantially high number of ion- 
izing photons and hence the ionized regions around them are rela- 
tively small. However, miniquasars of higher masses like 10^ and 
lO*' solar masses produce ionized bubbles of a couple of comoving 
megaparsecs. 

An interesting detail apparent in these Figures is that the Hel I- 
front exceeds that of the HI. Explanation for this lies in the fact that 
we have relatively high energy photons (E > 50 eV) and the cross- 
section of HI is much smaller than Hel, cTHeiiE) /aHi{E) ~ 20 
at an energy E = 50 eV. Therefore, the probability of high energy 
photons being capture by Hel is higher, increasing the helium I- 
front distance from the centre. 

Results are shown for the HE {SJ and LE (O cases. As ex- 
pected the ionized fronts of HI and Hel have travelled a much 
greater distance in the LE case. This is consistent because the HE 
spectrum is not only devoid of hydrogen ionizing photons but also 
lacks helium ionizing photons (54.4 eV). 

Another point to notice, is the visible excess of HI fraction 
in the HE spectrum just before the full ionization front in Figure|4] 
especially at the 10 Myrs curves. We interpret this excess as a result 
of the interplay between the increase of Hell and decrease of Helll 
as seen in Figure|6]thus increasing the HI recombination rate on the 
one hand and decreasing photon flux as a function of radius on the 



Figure 5. Same as in figE] but now including UV ionizing photons 
(10.4eV < E < 10-*oV). 

other. This interpretation is supported by the weakening of the HI 
excess feature as a function of black hole mass. The phenomena is 
not manifested in the LE spectrum (FigurelS} due to the ionization 
efficiency of the UV photons. 

4.3 The kinetic temperature profile 

The heating and cooling terms included in equation [12] are coupled 
with the rate equations [T] to [3] The principal heating terms are the 
bound free and Compton heating. Bound free heating is by far the 
most dominant, although towards the centre Compton heating does 
become important. 

Figures [8] and |9l corresponding to the HE and LE cases re- 
spectively, show the temperature profile for masses from WAIq to 
IO^'A/q for four different redshifts. The snapshot is taken after 3 
Myrs of evolution. 

Although there is a substantial change in the ionization pro- 
files (sizes of the ionized spheres), the heating remains more or less 
the same except for a few details. The flux of high energy photons, 
those responsible for the dominant secondary heating is abundant 
causing a sort of invariance in the extent of heating in both cases. 
Miniqsos of high masses do maintain an extended volume of high 
temperatures, in some cases up to 5 comoving Mpcs. If this indeed 
is the case at relatively low redshifts like around six, then they will 
be observable in the Sunyaev-ZeTdovich effect as discussed in sec- 
tion f|7] 

The four panels corresponding to different redshifts do not 
show a significant difference in the heating or its extent. Whereas 
the ionized bubbles (Fig. |4] and/or |5j are strongly dependent on 
the redshift. This is due to the fact that lower energy photons (re- 
sponsible for ionizing) see a much higher optical depth in neutral 
hydrogen than do high energy photons. Thus an increase in density 
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Radial Distance (l<pc) 

Figure 6. Fraction of various species (xhi,3:^hii.^HcI:2;hcII. 3;'hcIII ) 
as function of distance from tlie centre of the black hole with masses 
]^q1,2,...6^^ (top-left to bottom right panel) at a redshift z = 20. These 
are snapshots after 10 Myr of evolution. Miniquasars have no UV ionization 
energy (200eV < E < lO^'cV). 

by an order of magnitude (from z=10 to z=25) does not alter the 
heating considerably. 



4.4 Influence of different spectral energy distributions 

Throughout our study in this paper we have dealt with ct = 1 case 
for equation[T8] Here we explore the possibility of other indices for 
the power law. As an example we plot (ref. Fig. llOt the size of the 
ionized bubbles after 10 Myr as a function of a around a quasar 
with central black hole of lO^i\f0 for four different redshifts as 
indicated in the figure. Fi sure [Tol shows that the size of the ionized 
bubbltO increases with a. This is expected because an increase in 
Q for the same given normalization concentrates more photons at 
the low energy end, i.e., around 13.6 eV and fewer high energy 



^ We define the size of the ionized bubble as the position of the I-front 
when the ionized fraction is 0.5 for hydrogen 



Radial Distance (l<pc) 

Figure 7. Same as in fig|6] but now including UV ionizing photons 
(10.4eV <E < lO-^eV). 

photons. Thus, for higher a, one expects more ionization and less 
heating. This indeed is reflected in Figure [TTI 

The central regions are heated more or less to the same extent. 
In order to see the significant difference we have only plotted the 
distances between 100 and 1000 kpc physical. It is clearly seen that 
for steeper slope heating drops considerably. 



4.5 Life-time & Duty-cycle of mini-quasars 

IWvithe & Loebl ( l2002h produce analytical estimates of the quasar- 
life time and duty cycle (see their Figure 3). Certainly this sort of 
periodic switching on and off of the quasar leaves an imprint on the 
IGM. We simulate a case wherein the quasar duty cycle is 100 Myr 
with an active phase of 10 Myr at the start of the cycle. 

Recombination time-scales are orders of magnitude larger 
than the ionization timescales at the mean densities of the IGM. As 
a result, once the IGM is ionized by a source with finite life time 
(say lOMyr), it takes considerably long time for the IGM to recom- 
bine, leaving behind a bubble of ionized gas. Figure [12] shows the 
ionized fraction and the temperature profile just after the source is 
switched off at 10 Myr, than at 50 Myr and at the end of the duty 
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Figure 8. Radial profile of Kinetic temperature Te for black hole with 
masses 10'^'^'--^Mq at 4 different redshift z = 10, 15, 20, 25 (top-left 
to bottom right panel) after 3 Myrs is shown. The spectrum of these quasars 
include only high energies. 



cycle which is 100 Myr. The source is a miniqso of a 1000 Mq at 
redshift 20. 

There are a couple of points to be noted here. One is the fact 
that the I-front position is further away at 50 and 100 Myr than 
at 10 Myr although the source was switched off after 10 Myr. This 
propagation of the I-front in the absence of photons from the source 
is attributed to the ionization due to collisions. This continues until 
the temperature drops sufficiently so that collisions become inef- 
fective. In fact before the IGM recombines completely the quasar 
is switched on again. Thus we seeQthat the ionizing front expands 
during the on-time, stands still around the same position during the 
off-time and continues to expand during the next on-time and so 
on. 



^ movies created from the simulation can be obtained by contacting the 
authors 
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Figure 9. Same as in fig[8]but these quasars do include the low energy 
photon. 

The drop in temperature on the other hand is easily seen within 
this turn-off period. Compton cooling is the dominant sink for the 
temperature. Also we see that the regions close to the centre remain 
relatively more ionized than the rest. This is again a consequence of 
the temperature profile. Regions with higher temperature continues 
to be ionized for longer periods. 

All of the analysis performed thus far was done by embed- 
ding the source in a uniform background. But we know that objects 
like quasars are preferentially formed in overdense regions. In order 
to test the impact of different overdensities on the I-front position 
and temperature produced by a quasar we embedd the quasar in a 
background that is 1,10 and 100 times the mean IGM density. The 
results are shown in Figure [T3] The miniqso has a 1000 Mq black 
hole at the centre and is at a redshift of 10 and contains only the 
high energy photons. As expected the ionization fronts could pen- 
etrate less at higher densities. But the recombination rate is much 
greater at higher densities. This is seen in the top right panel cor- 
responding to 100 Myr as the I-front corresponding to a clumping 
factor of a 100 has recombined much more. On the other hand the 
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Figure 10. Position of the I-front is plotted in physical size as a function of 
spectral index a, 10 Myr after the quasar was switched on . Note that the 
sizes of the ionized region plateaus after about a = 2.5. The black hole in 
the centre is of 10^ Mq. 
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Figure 11. Temperature around a 10® Mq black hole is plotted for three dif- 
ferent spectral indices at z=20. We see the decrease in heating with steeper 
indices. Decrease of low energy photons being the cause of this behaviour. 
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Figure 12. Figure shows the neutral fraction (top panel) and kinetic temper- 
ature (bottom panel) for three different times after the quasar was switched 
on as indicated in the figure. The quasar is at redshift z = 10 and hosts a 
1000 Mq black hole in the centre, which is switched on for the first 10 Myr 
and switched off since then. 

heating is relatively less at higher densities because most of the 
photons are absorbed to ionize and the dominant cooling terms, 
e.g., like the Hubble cooling, are enhanced at higher densities. 



5 APPLICATION TO STELLAR (BLACKBODY) 
SOURCES 

Unlike quasars and miniqsos that have a power-law distribution of 
energy in their spectrum, the stars approximately behave as a black- 
body of a given temperature. This blackbody nature of the source 
leaves different signatures in the manner the IGM is heated and ion- 
ized. Give n the wea k dependence of the star's temperature on its 
mass (see lSchaerejJIoOZ) ), we fix the blackbody temperature of 
the stars to 5 x 10* K and mass ranges between 10 and 1000 Mq. 
The tota l luminosity for a given stellar mass is calculated from ta- 
ble 3 of lSchaerej ( |2002|) . The other difference in case of stars is 
the density profile in which it is embedded. We assumed a density 
profile of the form; 

p{r) = 3.2 X (91.5 pc/rf [cm~^] (r < r^) 

= n,(0) X (1 + zf [cm~'^] {r > r^), (21) 

where 71^(0), is the density of the hydrogen or helium at redshift 
and rc is the radius at which the density profile falls to that of the 
mean IGM. 

Depending on the temperature, the spectrum used above for 
the stars peaks in between « 20 eV to « 24 eV (for 10^ K). What 
is expected in the case of stars of the same power output as a quasar 
are bigger ionized bubbles and less heating. Simply because, after 
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Figure 13. Figure shows the neutral fraction (top panels) and kinetic tem- 
perature (bottom panels) for three different clumping factors at 10 Myr and 
100 Myr as indicated. The quasar is active at redshift z = 20 and hosts a 
1000 Mq black hole in the centre, which is switched on for the first 10 Myr 
and switched off since then. 



Figure 14. The neutral hydrogen fraction as a function of distance from 
the centre of the blackbody with an effective temperature of 50000K and 
masses 10, 100, 500 , 1000 A/q (top-left to bottom right panel) at a redshift 
z = 20. The three lines corresponds to the I-front position after 2,3 and 5 
Myr of evolution. 



the peak, which is around the ionization threshold of hydrogen and 
helium, we have an exponential cutoff towards higher frequencies. 

After modifying the spectrum and the underlying density pro- 
file, we basically did the same exercise as for the miniqsos. Some of 
the results are plotted in the following Figures. The hydrogen neu- 
tral fraction is plotted in figure [14] as a function of distance from 
the source with an effective blackbody temperature of 50000 K. 
Note that unlike the miniqsos case the x-axis here runs only up to 
100 kpc physical. This is because the most massive stars consid- 
ered here, i.e., 1000 Mq, has a lumi nosity of a 10 Mq black hole 
shining at the Eddington rate jSchaerer 2002). 

A similar case of lower ionized helium is seen in Figure [Tsl 
Observe again that the helium front is leading. Because for high 
temperature stars the blackbody peaks at higher energies closer to 
the ionization threshold of helium. 

Figure [16] shows kinetic temperature as a function of radial 
distance for Population III stars with masses of 10 — 1000 A/q. 
Notice that even though the maximum temperatures reached are 
comparable to that of the miniqso, the extent and shape of the pro- 
file are distinctly different. There is a much sharper edge to the 
heating attributed to the fact the there are not many high energy 
(100 eV and above) photons and hence the mean free path of most 
photons is lower. 



6 OBSERVATIONAL EFFECTS 

There are many indications of an epoch of reionization in the Uni- 
verse. But observations currently are only able to provide us with 
either an integral limit (CMB data) or a lower limit (Gunn-Peterson 



troughs) on the redshift of reionization. Direct detection and study 
of this epoch currently rests on the future 21 cm radio observations 
and temperature fluctuations because of the SZ effect (detectable by 
planckEI) and probable direct observations of very high redshift 
Population III stars by JWST0. 

In this paper we will concentrate on the feasibility of radio 
observations on the detection and mapping of the epoch of reion- 
ization. A point to note is that from here on in the case of miniqsos 
we only use the low energy case because the heating would be less 
exaggerated and probably closer to reality. 



6.1 Spin Temperature 

The spin temperature Ts couples to the either the CMB temperature 
TcMB o r to the kinetic temperature Te, in the absence of othe r radio 
sources ( iFieldll 19581; iKuhlen. Madau. & Montgomervll2006h . Thus 
Ts can be written as a weighted sum between Te and Tcmb , 

rp ^ T* + Tcmb + VcoiTe + VcTe ^^2) 

1 + Vcol + IJa 

T* = hh'2icm/k = 0.0681 K, j/coi and j/q determines the effi- 
ciency of collisional and Lya coupling and is given by, 

Vcol = {C„ + a + Cp) (23) 



* http://www.rssd. esa.int/index.php?project=PLANCK 
" http://www.jwst.nasa.gov/ 



12 Rajat M. Thomas and Saleem Zaroubi 



1 .0 



0.8 



X 


0.6 


X 




X 




"a) 


0.4 


X 




X 




zn 
X 


0.2 


X 




X 






0.0 


c 


1 .0 


o 








o 




o 


0.8 


L_ 




1 

(A 




03 


0.6 


O 




<D 




Q_ 




(/! 


0.4 



0.2 - 



0.0 



1 M , 



100 M 




500 M . 



] OOOM 



: 



Radial Distance (kpc) 



Figure 15. Fraction of various species {xiij,xiiii,xyicI,xhcU, ^ihcIII ) tor 
tlie sliown star masses. 



s ^ is the spontaneous emis- 



where, ^lo = 2.85 x 10 
sion rate or the Einstein A-coefficient and Ch, Ce, and Cp 
are the de-excitation rates of the triplet due to collisions with 
neutral atoms, electrons, and protons, respectively. The empir- 
ical fits for these coefficients are identical to those used by 
iKuhlen. Madau. & M ontgomerv (2 00^, which is a combination o f 
results publ ished in Zvgelman OOO^lAllison & Dalgarnd jl969l) . 
iLisztl j200lh and [smith C1966) . And the Lya coupling coefficient 
(discussed in section 



Va = 



27AioTemeC 



(24) 



Here, Jo is the Lya flux density. For the miniqsos high energy 
photons Lya coupling is m ainly caused by collisional excitatio n 
due to secondary electrons (Chuzh ov. Alvarez. & Shapiro! |2003) . 
This process is accounted for by the following integral. 



Mr) = 



UHiir) / a{E)N{E-r)dE, 



(25) 



where fi2 — 0.416 is the oscillator strength of the Lya transi- 
tion, e & me are the electron's charge and mass, respectively. (I>a 
is th e fraction o f the absorbed photon energy that goes into excita- 
tion jShull & va n Steenberg 1985). The contribution of this term is 
important close to the miniquasar. 

Note that in reality, we have to include the frequency depen- 
dence of the Lya cros s section, the background cont inuum Lya 
photons and scattering l lChen & Miralda-Escudell2006l) . But since 
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Radiol Distance (kpc) 

Figure 16. Radial profile of Kinetic temperature Te for stars with masses 
10, 100, 500 and lOOOM© at 4 different redshift z = 10, 15, 20, 25 (top- 
left to bottom right panel) are shown. 



we are dealing only with "injected" photons at the Lya frequency 
we only include the cross section at the line centre with thermal 
broadening. 

Results of the spin temperature for a black hole of mass 
1O*M0 is shown in Figure [TT] for 4 different redshifts. For rea- 
sons discussed in the next section the spin temperature follows the 
kinetic temperature for considerable distance away from the source 
and then follows the background CMB temperature at larger dis- 
tances. 

The spin temperature for the case of a lOOOAf© at 4 different 
redshifts is shown in Figure (TS) The spin temperature, for the case 
of stars, is assumed to follow the kinetic temperature for a distance 
of 'ct', although this might be an over estimate because there is a 
sharp cut-off towards red ward of the black body peak reducing the 
number of Lya produced by the source itself. But with a combi- 
nation of the background Lya , it might be possible to couple the 
spin temperature to the kinetic temperature. 
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Figure 17. Radial profile of spin temperature Ts for a black hole with mass 
10'^ Mq at 4 different redshift z = 10, 15, 20, 25 (top-left to bottom right 
panel) for 1,2 and 3 Myrs of evolution after the source was switched on. 
Note that the curves identical to that of the kinetic temperature in Fig|9]ex- 
cept the spin temperature climbs back to Tcmb after about 1 000 kpc corre- 
sponding to the photon propagation distance in 3 Myrs. These miniquasars 
include the low energy photons. 



6.2 Lyman a coupling 

Lya coupling is known to be the most efficient mechanism to 
couple the kinetic temperature to the spin temperature, at least in 
the mean IGM. Collisional coupling only becomes important in 
overdense regions with gas ove r-densities 5 > 20[(1 + 2)/10]~^ 
jlUev et alj2003l : lFurlanetto & Loeb 2002). Since we are interested 
in studying brightness temperature fluctuation relatively far away 
from the source, the dominant mode of coupling for our purposes 
is Lya pumping. 

There are a number of ways to generate Lya photons. 1) The 
source itself produces continuum photons between Lya and Ly/3, 
which is then redshifted at a distance from the source to Lya. 2) 
Photons above the Ly/3 produced by the source pumps the elec- 
trons to higher levels which then cascade back to produce Lya. 3) 
Secondary electrons can deposit a part of their energy toeards excit- 
ing the hydrogen atom (Shull & van Steenberg 1985). 4) If we are 
looking at sources at lower redsh ifts (< 20) then there is enough 
background Lya from "first stars" ICiardi & Madaul ( l2003h . 

The Lya flux density estimated in equation |25] is efficient 
to couple the spin temperature to the kinetic temperature. The 



Figure 18. Radial profile of spin temperature Ts for stars with a mass 
lOOOAf© at 4 different redshifts z = 10, 15.0, 20, 25 (top-left to bot- 
tom right panel) are shown. The temperature of the star is assumed to be 
50000K. 



Wouthuysen-Field effect | | is efficient within the photon propa- 
gation (or light travel) time if these photons come directly from the 
source within the energy range of ionization (13.6 eV) and Lya 
(10.4 eV). But, even in the case when the source does not produce 
Lya photons, the secondary excitations are efficient enough to cou- 
ple Ts with Te to a large distance away from the source as seen from 
Figure [771 But at lower redshifts the background Lya flux would be 
large enou gh to de-couple Ts from Tcmb much further away from 
the source i lCiardi & Madai3l2003h . In the case of stars we do not 
expect a major influence of the secondary electrons on ionizations 
or excitations. Thus we assume that the Lya flux originates from 
the source/background. Hence, we allow the spin temperature to 
follow the kinetic temperature to the light travel distance "ct". 



6.3 Brightness temperature 

Under the "low-frequency" or Rayleigh- Jeans' approximation, the 
equation for blackbody radiation reduces to, 



i{v) = 



(26) 



The Wouthuysen-Field mechanism basically is the mixing of the hyper- 
fine state in neutral hydrog en in its ground state via interm ediate transitions 
to the 2p and above states jWouthuvsenll952tlKei3ll958h 
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Where is the intensity of radiation at frequency i/, k the Boltz- 
mann constant, c the speed of light and Tb is the brightness temper- 
ature. And in radio astronomy this limit is applicable. Brightness 
temperature can be measured differen tially as a deviation from the 
back ground CMB temperature Tcmb ( lFieldll95^ : ICiardi & Madaul 
120031) as. 



sn = 



26771^- (1 + S) xm ( 

( 



' l + z 
10 



0.3 



1 - 

1/2 



0.02 



(27) 



Where xm is the neutral hydrogen fraction (see Fig|4l(, 5 is the 
overdensity of hydrogen (atoms and ions) and Ts is the spin tem- 
perature (Ref Fig 1171 ). So given the cosmology and the redshift, 
which translates to a frequency of observation, the brightness tem- 
perature in equation |27]basically reflects the source characteristics 
in two ways. One that enters through the neutral fraction .thi, de- 
fined by the ionizing capabilities of the source and another through 
the spin temperature Ts which in effect reflects kinetic temperature 
dictated by the heating capacity of the source. Thus, we hope that a 
3D tomography could not only reveal a statistical signal, the power 
spectrum and the large scale structure of the epoch reionization but 
also give us clues on the nature of these first sources. 

As an example, we can consider the HE and the LE case. Al- 
though the heating pattern and extent in both these cases are similar 
(see Figure[8]&|9j, the ionization profiles are indeed different (see 
Figure|4]&|5j. This according to equation|27]will be reflected in the 
brightess temperature. 

Figure [T9l shows the radial profile for the expected brightness 
temperature for the same case as discussed in the context of the 
spin temperature above. For all cases there is almost a sudden jump 
in the brightness temperature. This is due to the sharp transition in 
the neutral fraction (Ref FigQ. 

The ionized fraction (Fig|14t and spin temperature profile 
(FigllSt is also reflected in the brightness temperature in the case of 
stars as in Figure |20l The maximum of the brightness temperature 
has a typical value between 20 and 30 mK for all redshifts but the 
negative wing spans different ranges depending on the redshift. 

An interesting observation is that the maximum value of Tt 
remains the same throughout the mass range. The reason being that 
in all cases Ta is driven far away from Tcmb, such that the term 
(^1 — 2£MB.^ in equation l27l approximates to unity thus removing 
the dependence of Tt on Ts. As a consequence of this, we also 
see that the brightness temperature does not change considerably 
during the dorment period of the quasar. Figure [T3l shows ionized 
fraction and temperature profile after the quasar has been switched 
off. But the ionized fraction is still low (~ 10^^) in the inner re- 
gions of the I-front, and the temperature is considerably high. Thus, 
the brightness temperature profile essentially remains the same. 

IZaroubi & Silkl Jiool) have found that for power-law sources 
like miniqsos the brightness temperature has an extended feature 
in the transition from the ionized to neutral IGM. This feature has 
the same order of magnitude of the one found here (see Figure[79ll. 
However, it should be noted that the ori gin of this feature in both 
cases is completely different. While in IZaroubi & Silkl ( |2005|) the 
extended feature is produced by the ionization-recombination bal- 
ance of high energy photons, here the origin of this feature is in 
the heating not ionization pattern around these sources, also caused 
by high energy photons. As Figure |4] shows the neutral fraction of 
hydrogen changes abruptly at the I-front location, namely, we here 
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Figure 19. Radial profile of brightness temperature dTb for the same black 
hole model as in fig |17l 



show that the equilibrium assumption of IZaroubi & Silkl ( |2005|) is 
never attained during the miniqso lifetime. 



6.4 Observability with future telescopes 

As mentioned before, there are several interferometric radio tele- 
scopes underway to measure the 21cm signature of the EoR. It 
would be of interest to investigate into the possiblity measuring di- 
rectly the topology of the 21cm through the brightness temperature. 
This of course depends on the angular resolution, that measures the 
topology on the sky and the frequency resolution that will provide 
the depth information. We look into these aspects in the context 
of LOFAR and SKA measuring the 21cm signal from redshift 10. 

will have an angular resolution of about 3-5 arcmin and 
a spectral resolution of IMHz. This would correspond to a phys- 
ical size on the sky of ~ 1.4 Mpc and ~ 5 Mpc in depth. On 
the other hand SKA with an angular resolution of about I arcmin 
jCarilli & Rawling3l2004l) would correspond to ~ .5 Mpc on the 
sky. The spectral resolution for the EoR experiment would more or 
less remain the same to beat other noise effects. In the light of these 
rough estimates we realize that it would be unlikely for LOFAR to 
see individual sources, although once the "bubbles" start merging 
the sizes may become comparable. But SKA, in principle should 
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Figure 20. Brightness temperature calculated for the spin temperatures of 
the star shown in fig llSI Assumption has been made that the spin tempera- 
ture remains coupled to the kinetic temperature for a radius "ct". 



be able to resolve at least the relatively big bubbles of the order of 
100 of kpc easily. 



7 SUNYAEV-ZE'LDOVICH EFFECT 

Possibilities of observing the thermal and kinetic Sunyaev- 
Ze'ldovich (SZ) eff ect during the ep och of reionization had been 
first invest i gated by lAghanim et alj (1996) and more recently by 
llliev et alj feoOTh . A hot ball of ionized gas around first stars or 
miniqsos resembles the classic case of clusters at lower redshifts 
in which the SZ effect is observed today. Thus looking for signs 
of reionization through the SZ effect is an obvious next step. A 
spin-off of the simulations of the ionization and heating of the IGM 
by miniqsos is the estimation of their contribution to the Sunyaev- 
Ze'ldovich (SZ) effect. Both, the thermal and the kinetic SZ-effects 
have been considered. 

For these simulations we start at z of 20 with two seed black 
hole masses, 10^ and 10*Afo. These two cases are explored in or- 
der to reach, through accretion, masses of quasars observed at red- 
shift 6 by the SDSS. The black holes are allowed to grow in time 
according the following equation. 



M{t) = Mo exp{fdutyt/tE), 



(28) 
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Figure 21. Thermal Sunyaev-Ze'ldovich effect: Radial profile of tempera- 
ture fluctuation for seed black hole masses 10^ Mq and 10* Mq be- 
ginning at redshift 2 = 20 are shown. The different lines corresponds to 
particular redshifts as indicated in figure. Corresponding mass according to 
eauation l28l is shown. 



is around 41 x (erad /0-1) Myrs. The radiation efficiency of e^ad 
was set to 0.05 l lLoebl i2006). We assumed a duty cycle of 30%. 

First objects are formed in regions that are overdense with re- 
spect to the IGM. Thus, we ran the above simulations for a density 
that is ten times the mean IGM density at that epoch. 
Thermal SZ-ejfect: 

Hot electrons transfer energy to the CMB photons via inverse 
Compton scattering, redistrib uting the photons in the spectrum gen - 
erally towards higher energy dSunvaev & Ze'ldovicri[l972lll980l) . 
The temperature fluctuation is thus given by; 



AT ^ ^ fcfcCTT 



Te{t;l)ne{t-l)dl 



(29) 



where M{t) is the mass of the black hole at time 't\ fduty the 
duty cycle of the black hole and ts the time scale of growth which 



where Te,me and Ue are the electron's temperature, rest mass and 
density, respectively, ctt is the Thomson cross section and R is the 
radius of the ionized region. Te and rie are functions of the radius 
and time after the quasar was switched on. All parameters required 
to compute equation l |29t have been obtained from the simulation. 

The SZ-effect is the line of sight integral of the electron pres- 
sure. In order to estimate the effect we considered the profiles of 
the temperature and the ionized fraction to be spherically symmet- 
ric and then projected the effect onto the plane of the sky. Figure 
[2T]shows the temperature fluctuations resulting from the black hole 
model mentioned above. Results are plotted for different redshifts 
along with their corresponding masses, down to a redshift of 6. All 
of them end at redshift 6 with masses comparable to the black hole 
masses estimated for the SDSS high redshift quasars jpan et all 
2003, 2006). Massive miniqsos can cause temperature fluctuations 
at high redshifts on the order of lO^^K. 
Kinetic SZ-ejfect 

Apart from the thermal SZ-effect, there is also the possibility 
of observing these large ionized bubbles through the kinetic SZ- 
effect. Peculiar Motion of an ionized bubble with respect to the 
background leads to an additional change of radiation temperature 
in its direct ion, because of the finite opti cal depth associated with 
the bubble l lSunvaev & ZeTdovich|[T973) . The average anisotropy 
associated with an ionized bubble is given by, 

^ = 1^ / 2aTn,[t-l)dl (30) 

where Vr is the radial component of the peculiar velocity of the ion- 
ized bubble. The electron density rie (t; I) is given by the simulation 
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and the radial component of the velocity dispersion is computed in 
the limit of linear theory and is given by, 

Vr(z) = VriO){l + zy^^\ (31) 

where Vr{0) is velocity dispersion at redshift zero fixed at 600 
km/s. 

The temperature fluctuation as a result of the kinetic SZ effect 
(eg. 1301) is shown as a function of the radial distance from the source 
in figure [22I Again, as for the previous case, we assumed spherical 
symmetry. The kinetic SZ- effect is in fact larger in magnitude than 
thermal SZ-effect, making it easily visible on the temperature maps 
of PLANCK. Therefore, if the sizes and temperatures of these early 
ionized bubbles are of the same order of magnitude predicted here, 
then the SZ-effect is an important and independent probe of the 
impact of early sources on the reionization history that could be 
observed in the near future. 



8 SUMMARY AND CONCLUSIONS 

In this paper a 1-D radiative transfer code has been developed in 
order to study in detail the influences of primordial black holes and 
Pop 111 stars on their surrounding environment. We have shown, in 
agreement with many other authors, that black holes and/or Pop 111 
stars have the potential to have been the primary source of ioniza- 
tion. Although the ionized regions around typical black holes and 
Pop 111 stars is not shown to be convincingly different modulo the 
resolution of the telescopes like LOFAR, the heating around these 
two sources do have a significantly different profile. This in turn 
is reflected in the spin temperature of the system, which then di- 
rectly translates into the brightness temperature measured by the 
radio telescope. 

Spin temperatures are coupled to the kinetic temperature ei- 
ther through collisions or through Lya pumping. Strengths of these 
coupling terms are studied as a function of the radial distance from 
the source. Lya photons do couple the spin temperature efficiently 
to the kinetic temperature hence making the brightness tempera- 
ture insensitive to the background CMB, at least for a substantial 
distance (~ 2 — 3Mpc) away from the source. 

Recombination timescales are orders of magnitude lower than 
that of ionization. Therefore, an Hll bubble remains ionized for 
a long period after the source of radiation has been switched off. 
Also the temperature remains largely unchanged for a significant 
fraction of Hubble time since cooling processes typically involve 



two body interaction that increase the timescales involved. Inter- 
estingly for miniqsos, a marginal increase in the ionized bubble 
after the switching off of the central ionizing source is observed. 
Collisional ionizations are the dominant cause of this phenomenon 
which is therefore relatively less in stars because of the lack of high 
energy X-ray photons which heats up the environment significantly 
hence boosting the collisional ionizations. 

For miniqsos with a hard spectrum that lacks the UV part of 
the spectrum we observe a excess of HI fraction as a function of 
radius just before the main ionization front. This is interpreted as 
a result of the interplay between the increase of Hell and decrease 
of Helll, which increase the HI recombination rate,on the one hand 
and the decrease in the photon flux as a function of radius on the 
other. This interpretation is supported by the weakening of the HI 
excess feature as a function of black hole mass. This phenomenon 
is not observed in miniqsos with spectrum that has UV photons, 
which is due to the ionization efficiency of these photons. 

The results of the IGM heating around miniquasars have als o 
been compared to the analytical approach of lZaroubi et alj ( l2007h . 
Results shown in Figure IaTI demonstrate that the analytical solu- 
tions, while in agreement to within an order of magnitude, under- 
estimate the kinetic temperature away from the centre by roughly a 
factor of up to 5. The difference in the results is due to the helium 
species cross-sections are not being fully accounted for in the an- 
alytical approach. In fact, in order for the heating to be accurately 
calculated one needs an exact knowledge of the abundance of each 
of the hydrogen and helium species averaged over their evolution 
history until the point in time in which one is interested and such 
detailed knowledge could not be obtained with the analytical ap- 
proach. 

Also, the thermal and kinetic Sunyaev-Ze'ldovich effects 
around these quasars were studied. Two cases are considered in 
which mass growth of the black hole is incorporated such that it 
winds up with 10* and 10^ M0 black holes at redshifts close to six. 
The estimated values for the temperature fluctuations are within the 
sensitivities of the future mission like PLANCK. 

We conclude that the brightness temperature which in turn 
reflects the underlying ionization and heating around a source is 
sensitive to many factors. The spectral energy distribution of the 
source, range of energies spanned by the photons of these sources, 
their clustering properties, the photon escape fraction, the redshifts 
at which the sources turn-on, their lifetimes, and many other com- 
plex feedback mechanisms, to name a few that will influence the 
brightness temperature, which is the observable we are after. Thus, 
in order to span this large parameter space of possibilities a ID 
radiative transfer code is very useful. 

An important next step for the radiative transfer code de- 
scribed here is to include line emissions and absorptions, especially 
that of the Lya line, in the radiative transfer equations. We will also 
incorporate the rate equation for molecular hydrogen and its vari- 
ous states. 

Finally, in the near future we plan to incorporate the results of 
this study into the output of N-body simulations. Such an approach 
will allow us to produce quick 21 cm maps of the EoR, which then 
will be used in conjunction with the simulations of galactic and 
extra-galactic foregrounds, ionospheric models and LOFAR spe- 
cific instrument responses to generate "dirty-maps" of the EoR. 



Time-evolution of ionization and heating around first stars and miniquasars 17 



ACKNOWLEDGEMENTS 

The authors thank B. Ciardi, A. Nusser , E. Ripamonti and M. 
Spaans for discussion and helpful comments. We are also thank- 
ful to the anonymous referee for his illustrative and constructive 
comments. 



REFERENCES 

Abel T, Wise J. H., Bryan G. L., 2007, ApJ, 659, L87 
Aghanim N., Desert F. X., Puget J. L., Gispert R., 1996, A&A, 
311, 1 

Alhson A. C, Dalgarno A., 1969, ApJ, 158, 423 
Barkana R., Loeb A., 2001, PhR, 349, 125 

Begelman M. C, Volonteri M., Rees M. J., 2006, MNRAS, 370, 
289 

Benson A. J., Sugiyama N., Nusser A., Lacey C. G., 2006, MN- 
RAS, 369, 1055 
Bharadwaj S., Ali S. S., 2005, MNRAS, 356, 1519 
Bromm V., Kudritzki R. P, Loeb A., 2001, ApJ, 552, 464 
Carilli C. L., Gnedin N., Furlanetto S., Owen P., 2004, NewAR, 

48, 1053 

CarilU C., Rawlings S., 2004, astro, |arXiv:astro-ph/0409"274 
Chen X., Miralda-Escude J., 2006, astro, "arXiv: astro-ph/0605439 1 
Chuzhoy L., Alvai'ez M. A., Shapiro P R., 2006, ApJ, 648, Li 
Ciardi B., Ferrara A., White S. D. M., 2003, MNRAS, 344, L7 
Ciardi B., Madau P, 2003, ApJ, 596, 1 
Dijkstra M., Haiman Z., Loeb A., 2004, ApJ, 613, 646 
Dijkstra M., Haiman Z., Rees M. J., Weinberg D. H., 2004, ApJ, 
601, 666 

DopitaM. A., Sutherland R. S., 2003, adu..book, 
Elvis M., et al., 1994, ApJS, 95, 1 
Fan X. et al., 2003, ApJ, 125, 1649 
Fan X., et al., 2006, AJ, 131, 1203 
Field G. B., 1958 Proc I.R.E., 46, 240 
Fukugita M. , Kawasaki M., 1994, MNRAS, 269, 563 
Furlanetto S. R., Loeb A., 2002, ApJ, 579, 1 
Furlanetto S. R., ZaldaiTiaga M., Hernquist L., 2004, ApJ, 613,16 
Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1 
Hogan C. J., Rees M. J., 1979, MNRAS, 188, 791 
Iliev L T, et al, 2006, MNRAS, 371, 1057 
Iliev I. T, Pen U.-L., Bond J. R., Mellema G., Shapiro P R., 2007, 
ApJ, 660, 933 

Iliev I. T., Scannapieco E., Martel H., Shapiro P R., 2003, MN- 
RAS, 341, 81 

Jones B.J., Wyse R.F.G., 1985, A&A, 149, 144 

latayama T, Yoshida N., Susa H., Umemura M., 2004, ApJ, 613, 
631 

Kuhlen M., Madau P, Montgomery R., 2006, ApJ, 637, LI 

Kuhlen M., Madau P, 2005, MNRAS, 363, 1069 

Laor A., Fiore F, Elvis M., Wilkes B. J., McDowell J. C, 1997, 

ApJ, 477, 93 
Liszt H., 2001, A&A, 371, 698 
Loeb A., 2006, astro, arXiv:astro-ph/0603360] 
Madau P, Efstathiou G., 1999, ApJ, 517, L9 
Madau P, Meiksin A., Rees M. J., 1997, ApJ, 475, 429 
Mellema G., Iliev I. T., Alvarez M. A., Shapiro P R., 2006, NewA, 

11, 374 

Nusser A., 2005, MNRAS, 359, 183 

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P, 
1992, Numerical Recipres in C. The Art of Scientific Comupting, 
Cambridge University Press, Cambridge. 



Ricotti, M., Ostriker, J. P, 2004a, MNRAS, 350, 539. 
Ricotti, M., Ostriker, J. P, 2004b, MNRAS, 352, 547 
Rhook K. J., Haehnelt M. G., 2006, MNRAS, 373, 623 
Sazonov S. Y., Ostriker J. P, Sunyaev R. A., 2004, MNRAS, 347, 
144 

Schaerer D., 2002, A&A, 382, 28 

Scott D., Rees M. J., 1990, MNRAS, 247, 510 

Shapiro P R., Giroux M. L., 1987, ApJ, 321, L107 

Shaver P A., Windhorst R. A., Madau P, de Bruyn A. G., 1999, 

A&A, 345, 380 
Smith F J., 1966, P&SS, 14, 929 
Spaans M., Silk J., 2006, ApJ, 652, 902 
Spergel D. N., et al., 2007, ApJS, 170, 37 
ShuU J. M., van Steenberg M. E., 1985, ApJ, 298, 268 
Sunyaev R. A., Ze'ldovich I. B., 1980, MNRAS, 190, 413 
Sunyaev R. A., Zeldovich I. B., 1975, MNRAS, 171, 375 
Sunyaev R. A., Ze'ldovich Y. B., 1972, CoASP, 4, 173 
Tozzi P, Madau P, Meiksin A., Rees M. J., 2000, ApJ, 528, 597 
Vanden Berk D. E., et al., 2001, AJ, 122, 549 
Vemer, D. A., Ferland, G. J., Korista, K. T, & Yakovlev, D. G. 

1996, ApJ, 465, 487 
Vignali C, Brandt W. N., Schneider D. P, 2003, AJ, 125, 433 
Volonteri M., Lodato G., Natarajan P., 2007, arXiv, 709, 

arXiv:0709.0529 
Whalen D., Abel T., Norman M. L., 2004, ApJ, 610, 14 
Wouthuysen S. A., 1952, AJ, 57, 31 
Wyithe J. S. B., Loeb A., 2004, ApJ, 610, 1 17 
Wyithe S., Loeb A., 2003, astro, arXiv: astro-ph/ 0302297 1 
Wyithe J. S. B., Loeb A., 2002, ApJ, 581, 886 
Zaldarriaga M., Furlanetto S. R., Hernquist L., 2004, ApJ, 608, 

622 

Zaroubi S., Thomas R. M., Sugiyama N., Silk J., 2007, MNRAS, 
375, 1269 

Zaroubi S., Silk J., 2005, MNRAS, 360, L64 
Ze'ldovich Y. B., Sunyaev R. A., 1969, Ap&SS, 4, 301 
Zygelman B., 2005, ApJ, 622, 1356 



APPENDIX A: ANALYTICAL APPROXIMATION 

A number of papers have outlined analytical approximations for 
the ionized fraction a nd heating of the IGM by min iquasars (cf., 
iMadau. Meiksin. & R ees ( 1997): IZaroubi et aljHoOTh ). In order to 
test the validity of these analytical solutions, we compare such a 
solution with the RT simulations for 100 and 10* M0 black hole 
masses at z of 20 and 10 and with a lifetime of 3 and 10 M yr. Here 
we follow the analytical equations of lZaroubi et alj l l2007h who cal- 
culate the ionization species by solving the following ionization- 
recombination equilibrium equation; 

alrjUuil - XHi) ^T(E;r)nHXHi . (Al) 

V (JH riH J 

where T{E;r) is the ionization rate per hydrogen atom for a given 
photon energy at distance r from the source. F is calculated as a 
function of r; 

T{E-r)= / a(E)M{E-r)'i^. (A2) 

a\jj is the recombination cross-section to the second excited 
atomic level and has the value 2.6 x 10""r4~° *''cm^s"\ with 
Ti, the gas temperature in units of 10* K,. For this calculation they 
assume T = lO* K. Although not very accurate, it gives a lower 
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limit on the recombination cross-section, Oj^'j (in neutral regions 
atomic cooling prevents the gas from having a higher temperature). 

Similarly the heating rate Ti.{r), is calculated as a function of 
the radial distance r using; 

n{r)=fnHXHi a{E)Af{E;r)dE (A3) 
Jeo 

where / is the fraction of the absorbed photon energy that goes 
into h eating through collis ional excitations of the surrounding ma- 
terial dShull & van Steenbe rg 1985). The function / is fitted in the 
paper by IShuU & van Steenberg! ( I1985t) with the following simple 
fitting formula: / = C [l - (1 - x")''], where C = 0.9771, 
a = 0.2663 and b = 1.3163 and x — 1 — xui is the ionized 
fraction. This fitting function is valid in the limit of high photon 
energies, an appropriate assumption for the case at hand. The fit- 
ting formula is modified by imposing a lower limit of 11% for the 
fraction of energy that goes into heating as the proposed fitting for- 
mula does not work well at ionized hydrogen fractions smaller that 
10"''. And a{E) is the cross-section combining that of hydrogen 
and helium in the following manner; 



a{E) = aniE) + !^aHe{E) = ai(^) 

UH \ E J 



(A4) 



Here, cri is a smooth function of energy. The temperature of the 
IGM due to this heating is determined by the follo wing equation 
dMadau. Meiksin. & Reeslll99l IZaroubi et al]|2007h : 



3 nnkbTkinir) 



= n{r) X 



(A5) 



Here T^in is the gas temperature due to heating by collisional pro- 
cesses, fci, is the Boltzmann constant, ^ is the mean molecular 
weight and tq is the miniquasar lifetime. This equation assumes 
that the heating rate due to the absorption of X-ray photons during 
the miniquasar lifetime is constant. Given the miniquasars lifetime 
of 3 and 10 Myr and Hubble time at the redshifts we are interested 
in, cooling due to the expansion of the Universe can be safely ne- 
glected. A lO^'K cutoff on the gas kinetic temperature due to atomic 
cooling is imposed. Figure IaTI shows the kinetic temperature pro- 
files for the 4 cases namely, redshifts of 20 and 10 and lifetimes of 
3 and 10 Myrs. 

Although the figure shows an overall agreement between the 
simplistic analytical model and the full RT numerical solution, 
the analytical approach underestimates the heating by a factor of 
roughly 3 to 5 in the outer radii. The reason for this is that the he- 
lium cross-section is accounted for in the analytical approach (see 
eq. IA4b in a very incomplete fashion whereby the heating stops as 
soon as hydrogen is completely ionized (eq. IA3l l. In reality this is 
not the case as heating by high energy photons will continue un- 
til also He is fully ionized. Figure lA2l shows the cross-sections of 
neutral hydrogen (dotted line), completely neutral helium (dashed 
line), singly ionized helium (dot dashed line) all weighted by their 
abundances. Over plotted on this is the case in which all hydrogen 
and half of all helium has been completely ionized (50% Hell case). 
At high energies, like the spectrum of the miniqso under discussion, 
even this fraction of singly ionized helium provides a substantial 
cross-section that can capture the photon and convert its energy to 
heat. Interestingly, looking at Figures [6] and |7] in conjunction with 
Figures [8] and [9] around a radial distance of 20 - 30 kpc reveals that 
there is a correlation between the increase in temperature and the 
corresponding increase in the singly ionized heliimi species. 
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Figure Al. Each panel in this Figure shows a comparison between the 
model of Zai'oubi et al. ( 2007) and the results of the spherically symmet- 
ric radiative transfer code applied to black hole masses, 100 & 10000 Mq 
with the same spectral energy distribution. The analytical calculation is rep- 
resented by the solid line and that obtained from the radiative transfer code 
is represented by the dashed line. The different panels coiresponds to dif- 
ferent values of redshift and quasar lifetime 'tq' 
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Figure A2. Figure shows the cross-sections <Ve meretal.]|l996l) of neu- 
tral hydrogen, neutral helium and singly ionized helium, all waited by their 
abundances in the mean IGM. Over plotted is the also the 50% singly ion- 
ized helium, i.e., all hydrogen has been ionized and 50% of all helium is 
completely ionized. 



